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Abstract 

We adapt and study a variance reduction approach for the homog¬ 
enization of elliptic equations in divergence form. The approach, bor¬ 
rowed from atomistic simulations and solid-state science [Ml [MIES], 
consists in selecting random realizations that best satisfy some statisti¬ 
cal properties (such as the volume fraction of each phase in a composite 
material) usually only obtained asymptotically. 

We study the approach theoretically in some simplified settings 
(one-dimensional setting, perturbative setting in higher dimensions), 
and numerically demonstrate its efficiency in more general cases. 

1 Introduction 

1.1 Overview 

In this article, we adapt, theoretically study and numerically test a 
specific variance reduction approach for the numerical homogenization 
of an elliptic equation with heterogeneous random coefficients. 

The equation we consider is the following scalar elliptic equation 
in divergence form 

— div Q Vri^(-,a;)^ = / in T>, tt^(-,a;) = 0 on PP, (1) 
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set on a bounded regular domain P in (for some d > 1), with a 
deterministic function / G in the right-hand side. The field 

A is a fixed matrix-valued random field. It is assumed to be uniformly 
elliptic, uniformly bounded and stationary in a discrete sense. All this 
is made precise in Section [1.21 Since the parameter e in ([1]) is assumed 
small, the coefficient A ^-, 0 ;^ is oscillatory and © is challenging to 
solve numerically. On the other hand, the problem is theoretically well 
understood, as is recalled below. 

In the numerical practice, the traditional approach to approximate 
the solution to © is to consider (for any p G M'^), and solve, 

the so-called corrector problem 

— div [A(p -|- VtCp)] = 0 in almost surely, 

I E(Vt(;p) = 0, VrCp is stationary in the sense of ([5]) below, 

. JQ 

( 2 ) 

associated to ©■ The solution to m gives the deterministic and 
constant coefficient A* of the homogenized equation that in turn serves 
for the approximation of ©■ We refer to Section [1.21 below for details. 

Since © is a problem set on the entire space M'^, it is necessary to 
truncate it on a bounded domain, and to complement it with appro¬ 
priate boundary conditions. In practice, it is standard to consider the 
problem 

-div ^A(-,a;) (p-l-Vt(;^(-,w)) ^ = 0, r(;^(-, w) is QAr-periodic, (3) 

where, say, Qn = The deterministic homogenized matrix A* 

is then approximated by the random variable A^(uj) defined by 

VpGM'^, Alf{uj)p=-^f A{-,u;){p + Vw^{-,uj)). (4) 

\Qn\ Jq^ 

This approximate homogenized coefficient A^(a;) is then evaluated 
using the Monte-Carlo method. Random realizations of the environ¬ 
ment, namely the matrix coefficient A(y,a;), are considered within the 
truncated domain Qn- For each of these environments, m is solved 
and the matrix A^{u}) is computed using ©. The homogenized coeffi¬ 
cient A* is eventually approximated as an empirical mean over several 
realizations of A^(uj). More details are given below in Section [1.31 

The purpose of this article is to reduce the variance of the approx¬ 
imation of A*. 
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For this purpose, we borrow a variance reduction approach orig¬ 
inally introduced in a completely different context, namely that of 
atomistic simulations for microscopic solid state science. In the se¬ 
ries of articles [231 (2^ I25| . an approach is indeed described that se¬ 
lects some particular random realizations of the environment, based on 
some selection criteria derived from asymptotic properties. Intuitively, 
the approach aims at considering only realizations that, for N fixed, 
already satisfy properties that are usually only obtained in the asymp¬ 
totic limit N —)• oo. The approach carries the name SQS, abbreviation 
of Special Quasirandom Structures. Its principles share some similarity 
with those underlying another classical variance reduction technique, 
namely stratified sampling. 

We aim at adapting this approach to our context, at studying it 
theoretically in some simple situations, and testing it numerically in 
more general situations. 

For the sake of completeness, we mention that we have already 
studied the theoretical properties and the practical performance of sev¬ 
eral variance reduction methods for numerical random homogenization 
in some previous works of ours. The classical approach of antithetic 
variables, an approach that is quite generic and does not require nor 
exploit knowledge of the specihc structure of the random problem at 
hand, has been considered in [I1151191II8]. The significantly more elab¬ 
orate (and thus more efficient) approach of control variates is the sub¬ 
ject of m That approach requires a better knowledge of the problem 
considered, and is not always amenable to fully generic situations. 

Our article is articulated as follows. 

In the remainder of this introductory section, we present the basics 
of the theoretical setting (in Section [LSI) and of the numerical approx¬ 
imation method (in Section fl.Sp for the homogenization of the random 
equation ([T]). 

In Section 121 we introduce the variance reduction approach we con¬ 
sider. For pedagogic purposes, we first briefly expose the approach in 
the context of solid state physics it has originally been introduced in. 
This is the purpose of Section 12.11 In Section 12.21 we formally derive 
the specihcs of our variance reduction approach using a perturbative 
setting. This formal derivation provides the motivation for the gen¬ 
eral so-called SQS conditions that we use in the sequel of the work. 
Section [2j3] presents how we compute these conditions in practice. Sec¬ 
tion 12.41 contains the pseudo-code of our approach, along with some 
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comments. 

The theoretical analysis of the approach is the purpose of Section [3j 
We begin by proving, in a fairly general situation (in any ambient 
dimension), that the approximation provided by our approach (at least 
the simplest variant of our approach) converges to the homogenized 
coefficient A* when the truncated domain converges to the whole space 
(see Theorem [8] in Section [3.ip . Next, in Section [3.21 we investigate 
more thoroughly the one-dimensional setting, where we can indeed 
completely analyze our approach and actually prove its efficiency. 

Our final Section 0] contains numerical tests. First, since it is of¬ 
ten necessary to enforce the desired conditions up to some tolerance 
(see Remark [3] below), we investigate in Section 14.11 how this toler¬ 
ance affects the quality of the approximation and the efficiency of the 
approach. We observe there that the approach is robust in this respect. 

In Section 14.21 we illustrate on a prototypical situation the effi¬ 
ciency of our approach and scrutinize its sensitivity and the various 
sources of error involved. The conclusions are the following. The sys¬ 
tematic error is kept approximately constant by the approach (it might 
even be reduced), while the variance is reduced by several orders of 
magnitude. The more conditions we impose on the microstructures, 
the smaller the variance. The total error is always reduced. Such an 
efficiency is achieved at almost no additional cost with respect to the 
classical Monte Carlo algorithm. 

In order to demonstrate the versatility of the approach, we ap¬ 
ply it in Section 14.31 to a case with a way more general geometry of 
microstructures. There again, the approach provides a significant re¬ 
duction of the variance. 

We conclude this overview by emphasizing that, although the ap¬ 
proach introduced in this article is applied to the simple linear elliptic 
equation there is no reason to believe that it cannot be applied for 
a large class of partial differential equations with random coefficients. 
Indeed, the principles of the approach do not depend upon the specihc 
form of the equation. 

1.2 Theoretical setting 

To begin with, we introduce the basic setting of stochastic homoge¬ 
nization. We refer to the seminal works dSlEI] , to m for a general, 
numerically oriented presentation and to [2113 IB] for classical text¬ 
books. We also refer to m and the review article [T] (and the extensive 
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bibliography therein) for a presentation of our particular setting. 

Throughout this article, (0, P) is a probability space and we de¬ 
note by E(X) = X{ui)dF{u!) the expectation of any random variable 

Jn 

X G L^{Q, dP). We next hx d G N* (the ambient physical dimension), 
and assume that the group (Z*^,+) acts on 0. We denote by 
this action, and assume that it preserves the measure P, that is, for 
all A: G Z'^ and all E & X, F{TkE) = P(ii'). We assume that the action 
r is ergodic, that is, if ill G is such that TkE = E for any k G Z'^, 
then P(ii') = 0 or 1. In addition, we dehne the following notion of 
stationarity (see |16|): a function F G (M'^, L^(n)) is stationary if 

VA: G Z'^, F{x + k,uj) = F{x,TkUj) a.e. in x and a.s. (5) 

In this setting, the ergodic theorem |22| can be stated as follows: 

Let F G (M*^, L^(n)) be a stationary random variable in the 

above sense. Fork = {ki,k 2 , ■ ■ ■ ,kd) G Z'^, we set |A:|oo = max \ki\. 

l<i<d 

Then 

(2n\iY ^ in almost surely. 

This implies (denoting by Q = (0,1)'^ the unit cube in that 
F (yj F{x,-)dx'^ in L°°(M^), almost surely. 

Besides technicalities, the purpose of the above setting is simply 
to formalize that, even though realizations may vary, the function F 
at point X G and the function F at point x + A:, A: G Z'^, share 
the same law. In the homogenization context, this means that the 
local, microscopic environment (encoded in the matrix held A in ([1])) is 
everywhere the same on average. From this, homogenized, macroscopic 
properties follow. 

We consider problem m, which we recall here for convenience: 

— div Vu'^(-,a;)^ = / in P, u^(-,a;)=0 on dV. 

The random matrix A is assumed stationary in the sense of ([5]). We 
also assume that A is bounded and coercive, that is, there exist two 
scalars 0 < c < (7 < oo such that, almost surely, 

||A(-,a;)||£,oo(]Rd) < C and V.^ G A{x,uj)(, > c a.e. 
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In this specific setting, the solution to © almost surely con¬ 

verges (when e goes to 0) to the solution u* to the homogenized prob¬ 
lem 

— div {A*Vu*) = f inV, n* = 0 on dV. (6) 

The convergence of to u* holds weakly in H^(T)) and strongly 

in L2(P). 

The homogenized matrix A* in ([6]) is deterministic, and given by 
an expectation of an integral involving the so-called corrector function, 
that solves a random auxiliary problem set on the entire space. It is 
given by 


Vp E M'^, 


A* p = E 



A{x,-) {p + VWp{x, ■)) dx 


(7) 


where we recall that Q = (0,1)'^ and where, for any vector p € M'^, the 
corrector Wp is the unique solution (up to the addition of a random 
constant) in L^(fl; Lj^q^(M'^)) with gradient in L^(fl; of the 

corrector problem ([2]). We have used the notation Ly^j£(]R'^) for the 
uniform space, that is the space of functions for which, say, the 
norm on a ball of unit size is bounded from above independently of 
the center of the ball. 


1.3 Numerical approximation of the homoge¬ 
nized matrix 

As briefly mentioned above, the corrector problem ([2]) is set on the 
entire space M'^, and is therefore challenging to solve. Approximations 
are in order. In practice, the deterministic matrix A* is approximated 
by the random matrix A^(u}) defined by (jT]), which is obtained by 
solving the corrector problem ([3]) on a truncated domain, say the cube 
Qn = (0,^^)*^. Although A* itself is a deterministic object, its prac¬ 
tical approximation A^ is random. It is only in the limit of infinitely 
large domains Qn that the deterministic value is attained. As shown 
in |6], we indeed have 

lim A^(a;) = A* almost surely. (8) 

N^oo 

As usual, the error A* — A^{oj) may be expanded as 

A* - A%{uj) = - E [A^]) + (e [A^] - A^(a;)^, (9) 
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that is the sum of a systematic error and of a statistical error (the first 
and second terms in the above right-hand side, respectively). 

A standard technique to compute an approximation of E [A^] is 
to consider M independent and identically distributed realizations of 
the held A, solve for each of them the corrector problem ([3]) (thereby 
obtaining i.i.d. realizations for 1 < m < M) and compute 

the Monte Carlo approximation 


E 


i^N)i 


tMC/ 

-'M 




M 


( 10 ) 


m=l 


for any 1 < i, j < d. In view of the Central Limit Theorem, we know 


that E 


(A^). ■ asymptotically lies within the conhdence interval 


- 1.96- 


1 /var 


i WSLT 


V 





Vm 


y/M 


with a probability equal to 95 %. 

For simplicity, and because this is overwhelmingly the case in the 
numerical practice, we have considered in ([3]) periodic boundary con¬ 
ditions. These are the conditions we adopt throughout our study. It 
is to be remarked, however, that other boundary conditions may be 
employed. Likewise, other slightly modihed forms of equation ([3]) may 
be considered. The specihc choice of approximation technique is mo¬ 
tivated by considerations about the decrease of the systematic error 
in ([9]). Several recent mathematical studies have clarified this issue. 
In addition, in the particular case of periodic boundary conditions ([3]), 
it has been recently established in mi Theorem 2] that the statistical 
error in ([9]) decays like while the systematic error in ([9]) scales 

as Both estimates have been established for the discrete 

variant of the problem. A similar decay of the statistical error has also 
been established for the continuous case we consider in the present 
article (see m Theorem 1] and [201 Theorem 1.3 and Proposition 
1.4]). 
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2 Variance reduction approach 

2.1 Original formulation of the SQS approach 

The variance reduction approach we elaborate upon in this article has 
been originally introduced for a slightly different purpose in atomistic 
solid-state science [2311241 [25] , 

In order to convey to the reader the intuition of the original ap¬ 
proach, we consider here a simple one-dimensional setting, which nev¬ 
ertheless illustrates the difficulties of a generic problem. We consider a 
linear chain of atomistic sites of two species A and B which interact by 
the interaction potentials Vaa, Vab and Vbb with obvious notation. 
For simplicity we consider only nearest neighbour interaction. The 
atomic sites are occupied by a single species randomly chosen between 
A and B. A typical random configuration of the “material” therefore 
reads as an infinite sequence of the type • • • ABBAAABBAAAA ■ ■ ■ 

In order to compute the energy per unit particle of that atomistic 
system, one has to consider all possible such infinite sequences, and 
for each of them its normalized energy 



( 11 ) 


where Xi denotes the species present at the i-th site for that particular 
configuration (W = A or B). The “energy” of the system is then de¬ 
fined as the expectation of (|lip over all possible configurations. Other 
quantities than m may be considered, or may be simultaneously 
considered. 

In practice, one considers a presumably extremely large, finite N, 
truncates the infinite sequence over the finite length 2N + 1, and com¬ 
pute 



2N + 1 


i=-N 


for many (say M, where M is also presumably large) configurations. 

The approach introduced in [231 (2^ [25] consists in selecting spe¬ 
cific configurations (Xj)_ 7 v<i<Ar of atomic sites that satisfy statistical 
properties usually obtained only in the limit of infinitely large N. 

The first such statistical property is the volume fraction, namely 
the proportion of species {A,B) present on average. If the sites are 
all occupied randomly with probability \/2 oi A and 1/2 of B (and 




assuming that all these random variables are independent), then obvi¬ 
ously the volume fraction of A is 1/2 and so is that of B. Then, one 
only consider truncated sequences (Xj)_jv<i<Ar that exactly reproduce 
that volume fraction. 

Similarly, again for such an evenly distributed proportion of A and 
B, the energy of the entire infinite system evidently reads as 

£ = ^ [Vaa + ‘^Vab + Vbb] 


(recall that we only consider nearest-neighbour interactions). Thus, 
one only considers truncated sequences (Xj)_Ar<j<jv which, in addi¬ 
tion to exhibiting the exact volume fraction, have an average energy 

Vxi_f_iXi which is equal to £. And so on and so forth for 

i=-N 


N 

2N + l'^ 


other quantities of interest. 


Mathematically, this selection of suitable configurations among all 
the possible configurations classically considered in a Monte-Carlo 
sample amounts to replacing the computation of an expectation by 
that of a conditional expectation. 

The simplistic model we have just considered for pedagogic pur¬ 
poses can of course be replaced by more elaborate models, with more 
sophisticated quantities to compute, and more demanding statistical 
quantities to condition the computations with. The bottom line of 
the approach remains the same, and we adapt it to design a variance 
reduction approach for numerical random homogenization. 

In the next section, we derive the appropriate conditions, which we 
call the SQS conditions, for our specific context. 


2.2 Formal derivation of the SQS conditions us¬ 
ing a perturbative setting 

The purpose of this Section is to formally derive the SQS conditions 
that we use in the sequel. Such conditions can be easily intuitively 
understood. We however believe it is interesting to (formally) derive 
them in a particular case. The case we proceed with is a perturbative 
setting (although, we emphasize it, the conditions will be employed in 
the full general, not necessarily perturbative, setting). 

We assume throughout this section that the matrix valued coeffi- 
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cient A in m reads as 


Arj{x,uj) = Co{x,uj) + rj x{x,uj)Ci{x,oj) 


( 12 ) 


for some presumably small scalar coefficient rj, where 

• Co and Ci are two stationary, uniformly bounded matrix fields, 

• C'o(-, w) — C'i(-, w) and C'o(-,ti;) + C'i(-,a;) are almost surely coer¬ 
cive, 

• X is a stationary scalar field with values in [—1,1], 

Under these assumptions, for any rj E (—1,1), the matrix A^j is sta¬ 
tionary, bounded and coercive. Intuitively, when rj is small, A^j is a 
perturbation of the matrix-valued field Co{x,oj). 

Remark 1. The expression m models e.g. a two-phase composite 
material, where the phases are modelled by the coefficients Cq and Ci, 
while X is the indicator function of the first phase. 

Let p E The corrector problem ([2]) reads, in this particular 
setting, as 

( - div UCq + r]xCi){p+ Vwr,)] = 0 in M'^, 

{ e/^ V., = 0, .s in the sense of ©, 

and the homogenized matrix ([7]) is given by 



(14) 


Q 


Note that, for the sake of clarity, we omit to write the dependency of 
Wri with respect to p. 

The truncated version of (I13p on the domain Qjv is 


f -div [(Co -h??xC'i)(P + Vrc^)] = 0 in Qn 
1 i':^) is QAT-periodic, 

and we approach the homogenized matrix (|14l) by 


(15) 



( 16 ) 
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2.2.1 Expansion in powers of rj 

As T] goes to 0, we may now expand {u}) and A* in powers of r]. 
This expansion is classical (see for instance mm)- We only provide it 
here for the sake of consistency. The corrector expands as 


Vwr] = Vwo + rjVui + r]'^Vu2 + o{r]‘^). (17) 

This expansion holds in L^(n; L^jjj£(R‘^)). The functions wq, ui and 
U2 appearing in the expansion are respectively defined by the following 
systems of equations: 


— div [C'o(p + Vrco)] = 0 in 

E / Vtco = 0, Vrco is stationary, 
JQ 


(18) 


and 


— div [CoVtii] 

E f Vui = 0, 

JQ 


div [xCi{p + Vwo)] 
Vwi is stationary, 


in 


- div [C0VU2] 

E [ Vu2 = 0, 

JQ 


div [yCiVui] in 
Vtf 2 is stationary. 


(19) 


Inserting the expansion ([12]) of and (|17l) of Wrj in (|14p . we obtain 

A* = A* + r/At + r?2A* + o(r?2), (20) 

with, for any p G 


AIp = E 
A*^p = E 
= E 


[ Cq{p + Vwo) 
Uq 


Uq 


xCiip + Vwo) 


+ E 


CoVui 


[ xCiVui 

Uq 


+ E 


UQ 

C'oVrt2 


( 21 ) 


'Q 


Likewise, we expand as 


= VtCg + pVui + r]vu2 + oiv ), 


with 

f -div [C'o(p +Vrc^)] = 0 in Qat, 
[ Wq (-, 0 ;) is Q^v-periodic, 


( 22 ) 
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f - div [CoVtif ] = div [xCi{p + )] in Qn, 

1 Ui{-,uj) is QAT-periodic, 

and 

f —div [C0VU2] = div in Qat, 

I U2 (^w) is QAf-periodic. 

The homogenized matrix [ui) therefore satisfies 


( 23 ) 




(w) - ^o’^(a;) + r]Al’'^ (w) + (‘^) < 


.*,N 


2 a*A / 


(24) 


where C is independent of r], N and w, and where the matrices ^g’ (cj), 
A^’'^(w) and A^^ { 00 ) are defined by 

[ Co{p + Vw^), 

\Qn\ Jq^ 

A];^{uj)p = [ xC'i(p + V'u;^) + i^ / C'oV?rf,(25) 

\Qn\ Jq^ Ivivi Jqjv 

[ xCiVui 

IVAf| JQjv 




CqVu^. 


2.2.2 SQS conditions 

In line with the motivation we have mentioned above in Section [H 
we are now in position to introduce the conditions that we use to select 
particular configurations of the environment within Qj\f for which we 
compute the solution to (IlSp . and, in turn, compute the approxima¬ 
tion (1161) of A*. Our conditions are based upon the comparison of (1211) 
and (I25p . 


Definition 2. For finite fixed N, we say that an environment cj € 14 
satisfies the SQS condition of 

• order 0 if Aq’^( co) = Aq, that is to say, for any p € 


IQivI Jqn 


Co{;u;){p+Vw^{;Uj))=E 


UQ 


Co(p + Vrco) 


(26) 


• order 1 if Af^{oj) = A\, 


that is to say, for any p G R'^, 


1 

\Qn\ 


f xi-,ACii-,oj){P + '^Wo (•,a;))-bC'o(-,a;)Vuf (-, 0 ;) 
d Qn 


= E 


/ xC'i(p +Vrco) TOoVui 

Jq 


, (27) 


12 



















• order 2 if {oj) = A 2 , that is to say, for any p G 



= E f xCiVui + CoVu2 ■ (28) 

Jq 


Remark 3. In full generality, we do not claim that there exist envi¬ 
ronments that satisfy these conditions. This might he the case that no 
such environment exists. One may for instance simply remark that a 
random variable that takes value —1 and +1 both with probability 1/2 
never has value zero, which is its expectation! In some situations, we 
therefore have to relax the above conditions (see Section lKf] below), but 
we temporarily leave these technicalities aside and assume that suitable 
environments exist. 

Consider now the two expansions (|20l) and (1241) . It is immediate to 
see, by subtraction, that 



Therefore it is readily seen that, if the configuration uj satisfies the SQS 
conditions of Definition [2] up to the order k included (A: = 0, 1, 2 in 
our definition, but clearly one could consider higher order conditions 
derived likewise), then 



(29) 


where the constant in the right-hand side is independent of r], N and 
UJ. Taking the expectation over such configurations therefore formally 
provides a more accurate approximation of A*. 

Now that we have derived the conditions (j26jl - (l27p - (l28jl (which we 
henceforth call the SQS conditions) in the perturbative setting, we will 
actually use them in the non perturbative setting, namely for a simi¬ 
lar two-phase composite material, but with rj not small. Of course, a 
property such as (j29p cannot be expected any longer since the homoge¬ 
nized matrix A* is no longer a series in a small coefficient that encodes 
a perturbation. Nevertheless, it can be expected that selecting the 
configurations using these conditions may improve the approximation, 
in particular by reducing the variance. We show in Sections [3] and [J] 
that it is indeed the case, theoretically and experimentally. 
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For the time being, we need to make a practical observation. The 
right-hand side of conditions (I26p - (j27l) - (l28l) need to be evaluated in 
order to practically encode the SQS conditions. In principle, those 
right-hand sides are exact expectations, that can only be determined 
using an asymptotic limit, and are therefore as challenging to compute 
in practice as A* itself. 

We therefore need to restrict the generality of our setting (I12p and 
consider cases where those right-hand sides are indeed amenable to 
a simple, inexpensive computation. This is the purpose of the next 
section. 

2.3 Practical evaluation of the SQS conditions 

In order to make our approach practical, we need, as mentioned above, 
to consider settings where the expectations present in the right-hand 
sides of (IMP-dITp-dIHI) may be computed effectively. 

2.3.1 Condition of order 0 

We first consider (j26p and its right-hand side 



(30) 


Q 


A natural assumption, which already covers a large portion of practi¬ 
cally relevant situations, is 


Co{x,U}) = Co(x) is a deterministic, Z'^-periodic matrix. (31) 


The computation of (j30p is then inexpensive since the solution wq 
to (1181) is in fact the deterministic solution to 

— div [Cq{p + Vtco)] = 0 in M'^, wq is Z'^-periodic, 

which is unique up to the addition of a constant. 

In addition, when N is an integer (and when the approximation 
chosen for ([2]) is the periodic approximation ([3]), as is indeed the case 
throughout this work), the solution to (f22P is Wq = rco (up to an addi¬ 
tive constant), and hence the condition (|26p is systematically satisfied. 

We henceforth assume that (j3ip holds, that N is an integer, and 
that we proceed with the periodic approximation (l3|). 
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2.3.2 Condition of order 1 

We next consider the SQS condition (|27p . One possible assumption to 
make that condition practical is 


Cq{x,u) = Cq is a deterministic, constant matrix. (32) 


Since Vtoo = 0, the right-hand side of (|27l) reads 


E 


,/ 

UQ 


xCi{p + Vmo) -t CoVni 


= [ E[xCi]p + CoE [ Vm, 

JQ Jq 


where the rightmost term vanishes in view of (I19p and where the first 
term of the right-hand side may be computed using only characteristic 
properties of the environment considered. The condition (|27p thus 
reads 


1 

IQa^I 



x(-,w)C'i(-,w) = E 



(33) 


For instance, in a two-phase composite material mixing two con¬ 
stant and deterministic matrices Co and Ci, we have 


E 




Cl. 


This quantity obviously only depends on the volume fraction of the 
two phases (recall (I12li ). Proceeding likewise with the left-hand side 
of the condition (|27p . we see that this condition reads 


1 

IQa^I 


jT 

'a) 

II 

f/^l 

' Qn 

IJQ \ 


Interestingly (and not unexpectedly), we notice here that this con¬ 
dition on the volume fraction agrees with the condition we used to 
consider in the simple atomistic system of Section 12.11 


2.3.3 Condition of order 2 

We next proceed with condition (I28p . In addition to (I32p . we assume 
that 

Ci(a:,a;) = Ci(x) is a deterministic, Z'^-periodic matrix, (34) 
and that 

x{y,^) = ^ Xk{uj)tQ+k{y), (35) 

feeZ'* 
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where are identically distributed scalar random variables taking 
their values in [—1,1]. We also assume that 


C= ^ |Cov(Xo,Xfc)| <oo, (36) 

which is obviously satisfied if Xf^ are independent one from each other. 

We then have the following result, which will be useful to make 
condition (I28p practical. Its proof is postponed until Appendix lAl 

Lemma 4. Under the assumptions (|32p . (IMP . ([MP and (I36p . the so¬ 
lution ui to dm) satisfies 

Vufiy,u) = E[Xo]Vfir(2/) + (Xfc(w) - E[Afc]) V(?ii(y - k), (37) 


where is the (unique up to the addition of a constant) solution in 
{v € Vv G to 

— div [CqX fii] = div [IgCip] in (38) 

and Ml is the (unique up to the addition of a constant) solution to 


— div [CoVmi] = div 


Cip 


in Ml is W/"-periodic. (39) 


The sum in dszp is a convergent series in LfiiQ x fl). 

Using simpler arguments, we see that the solution Ui to (I23P sat¬ 
isfies 

VMf(y,w) =E[Xo]VMr(y)+ (^Xk{u)-E[Xk])x4>f{y-k), 

fceZ'inQjv 

_ (40) 

where mi is defined by (I39p and is the (unique up to the addition 
of a constant) solution to 

— div = div [IqCip] in Qiv, is QAr-periodic. (41) 


In practice, we can easily obtain an accurate approximation of 
<fi since the right-hand side of (|38p has compact support. Truncat¬ 
ing ()38P over a sufficiently large bounded domain (with homogeneous 
Dirichlet boundary conditions) provides such an accurate approxi¬ 
mation. Given (j32P . the right-hand side of Condition (|28p rewrites 
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E 

[ xCiVui 

since E 

[ Vn2 


Uq j 


Jq 

is 

m turn expanded as 



0. In view of (1371) . this quantity 


E 


[ xCiVui 

Uq 

{E[Xo]f J CiVnl+^E J Xo(Xfc-E[Xfc])CiV0i(--fc) 


fceZ'^ 


C*iVui 


(E[^o])' [ 

JQ 

+ J] E J (Xo - E[Xo]) (Xfc - E[Xk]^CiVM- - k) 


(42) 


where, as mentioned above, V0i can be easily and accurately com¬ 
puted, while the series in A; G Z'^ may be truncated in an efficient 
manner because of the rapid decay at infinity of V(/>i (see [5l Lemma 
3.1]). 

We correspondingly expand the left-hand side of (1281) . The second 
term vanishes, while the first term reads, in view of (|40p . 


^ f x{y,^)Ci{y)Vu^{y,uj)dy 
>n\ Jqj^ 


+ 


IQivI 

E 

j&zdnQN 

E 

k,jez<inQN 


IQivI 


\Qn\ 


f Xj(a;)lQ+jC'iE[Xo]Vtti 
J Qn 

Xj{u)lQ+jCi(^Xk{u) - E[Xfc])V(/)f (• - k) 


C*i Vui 


+ 


+ 


+ 


(E[Xo ])2 f 

JQ 

E (VjH-e[x,]) ) a 
E iTWr / C,(XiM-E[Xj)v0f'(.-fc) 


CiXui 


E 


fceZ'^nQjv 

1 


k,j(ilAnQN 


\Qn\ 


'Q+j 


X,(a;) -E[X,])Ci(Xfc(a;) -E[Xfc])v<(- - k). 


(43) 
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In this particular (however still very generic) setting, we infer from (|42p 
and (|43p that Condition (I28p reads as 


\Qn\ 




fcjsQivnZ'i 


+ 


\Qn\ 


E[Xo] [Xkico)-K[Xk])l^ = Cov(Xo,X4/r, 

keg^nzd- 


kezd 


(44) 


where 


TOO _ 

-'fc — 

= 

^k,j 


jN _ 
-'fc — 


[ Ci{y)VMy), (45) 

■/ Q-\-k 

f Ci{y)X(pi{y - k)dy, (46) 

Q~^j 

[ Ciiy)XcP^{y-k)dy+ [ Ci{y)Xwiy)dy. (47) 
Jqm Jq 


2.3.4 Summary 


In the prototypical case where 

A{x,uj) = Co + xix,uj)Ci{x), 

where Co is constant, Ci is periodic and x takes the form (l35P (and 
where we consider the periodic approximation ([3]) of ([2])), we have 
that: 


• The condition (I26p (SQS condition of order 0) is systematically 
fulfilled. 


• In view of ([33]) and ()35p . the condition (I27p (SQS condition of 
order 1) rewrites as 


I 

\Qn\ 


Y, Xkiu;)=E[Xo]. 

fceZ'^nQiv 


(48) 


• In view of (|44l) . the condition (I28p (SQS condition of order 2) 
writes as 


I 

IQa^I 


kJsQffnZd 


+ 




■E[Xo] 




YCov{Xo,Xk)ir, (49) 

k£Zd 
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where Xi.{io) = Xj^{lo) — IE[Xfc]. 

The conditions (I48p and (149 p are henceforth called the SQS 1 and 
SQS 2 conditions, respectively. 


Remark 5. If (1481) is satisfied, then the coefficient in ()49P can be 
replaced by 



and there is no need to compute Ui. 

2.4 Selection Monte Carlo sampling 

We are now in position to describe the selection Monte Carlo sampling 
we employ. We recall that the classical Monte Carlo sampling reads 
as follows: 

Algorithm 1 (Classical Monte Carlo). 

For m = 1,..., M, 

1. Generate a random environment ujm- 

2. Solve the truncated corrector problem (j3p. 

3. Compute 



m=l 


In contrast, our selection Monte Carlo sampling algorithm, in the 
particular case described in Section 12.3.41 reads as follows: 

Algorithm 2. 

The algorithm requires a tolerance tol > 0, fixed by the user. 

1. Offline stage 

(a) Solve the equation (1381) . 

(b) Compute {I^)k£Z‘^ defined by (1451) . 

(c) Compute the right-hand side of the SQS conditions (I48p and (09]). 

(d) Solve the equations (I39p and (|4ip . 

(e) Compute defined by (06]) 

and (l47p . 

2. Online stage 

For m = 1,..., M, 
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(a) Generate a random environment Um- 

(b) Using and , compute the left-hand sides of (|48p and (j49l) . 

(c) If the left-hand sides differ from the right-hand sides by more 
than tol, return to Step \2a[ 

(d) Solve the truncated corrector problem (j3|). 

(e) Compute {tvm)- 


M 


Compute the approximation of A 


m=l 


Remark 6. As pointed out above, the series in k ^ in the right- 
hand side of (j4^ may be truncated in an efficient manner because of 
the rapid decay at infinity ofS/fi. Therefore only a few factors 
have to be computed at StevUR 


Remark 7. When several SQS conditions (in practice SQS 1 and 
SQS 2) have to be simultaneously satisfied, we simply add them up 
using some weighting parameter. We have not observed any particular 
sensitivity of our numerical results (collected in Section^ below) with 
respect to the adjustment of this parameter, provided it remains not too 
close to 0 and 1. 


We have already mentioned that, in many situations, there might 
not be any random environments that satisfy some, or all, of the SQS 
conditions (I26p - (j27p - (j28p we wish to enforce. Therefore, some adapta¬ 
tion is in order, and we have used in Algorithm [2] a tolerance parameter 
tol > 0 for the SQS conditions to be satisfied. 

However, if these conditions are enforced within some given toler¬ 
ance as in Algorithm [21 the following issue arises. Since the motivation 
for precisely considering the SQS conditions is that they are fulfilled 
asymptotically, the larger the truncated computational domain we con¬ 
sider (that is, the larger N), the less restrictive the conditions are, and 
therefore the less effective the variance reduction is likely to be. To 
circumvent this difficulty, a first possibility is to consider a tolerance 
that decreases when the size of Qn increases. We consider this vari¬ 
ant in our theoretical study of Section [3.2.11 below (see formula (|6^ ). 
More precisely, we require in Proposition 1141 that 

the SQS condition is satisfied with the tolerance — , 
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for some A. In practice, implementing such a threshold is not an easy 
matter, as the rate and the constants need to be adequately adjusted. 
In order to avoid such technicalities, we prefer to take a slightly differ¬ 
ent perspective, the purpose of which is to always select a fixed propor¬ 
tion of the original sample of the A4 environments drawn. Practically, 
we pick the M configurations that best satisfy the SQS conditions 
among the Ai configurations that have been drawn. 

The practical algorithm we employ is therefore as follows: 


Algorithm 3 (Selection Monte Carlo sampling). 

The algorithm requires a number of trials Ai, fixed by the user. 

1. Offline stage 1: same as the offline stage of Algorithmic 

2. Offline stage 2: selection step 

For m = 1,..., Ai, 

(a) Generate a random environment Um- 

(b) Using and , compute the left-hand sides of (|48p and (I49p . 

(c) Compute the error error^ between the left-hand sides and 
the right-hand sides of (I48p and (I49p . 

Sort the random environments (cOm)i<m<M according to error^. 
Keep the M best realizations, and reject the others. 


3. Online stage: resolution 

For m = 1,..., M, 

(a) Solve the truncated corrector problem ([3]). 

(b) Compute 


M 


Compute the approximation Tgqs ~ Jff cAff{ujm) of A 


m=l 


We wish to make a couple of comments about this selection Monte 
Carlo approach. 

In full generality, the cost of Monte Carlo approaches is usually 
dominated by the cost of draws, and therefore selection algorithms are 
targeted to reject as few draws as possible. 

In the present context, where boundary value problems such as ([3]) 
are to be solved repeatedly, the cost of draws for the environment 
is negligible in front of the cost of the solution procedure for such 
boundary value problems. Likewise, evaluating the quantities present 
in e.g. (I49p is not expensive. Therefore, the purpose of the selection 
mechanism is to limit the number of boundary value problems to be 
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solved, even though this comes at the (tiny) price of rejecting many 
environments. This also explains why we employ a simplistic rejection 
procedure for the selection, while in other situations of Monte Carlo 
samplings, one would invest in a more clever selection procedure. 

A second observation is that, as potentially for any selection pro¬ 
cedure, our selection introduces a bias (i.e. a modification of the sys¬ 
tematic error in Q). The point is to ensure that the gain in variance 
superseeds the bias introduced by the variance reduction approach. 

Our next section addresses some theoretical aspects of our ap¬ 
proach. 


3 Elements of theoretical analysis 

This section contains some elements of analysis that we are able to pro¬ 
vide. We begin with a (somewhat) general result of convergence, and 
next, in some simplified cases, study our approach more thoroughly. 


3.1 Proof of convergence of the approach 


Formally, our approach consists in replacing an empirical average pro¬ 
vided by the classical Monte Carlo approach to compute E[A^] by an 
empirical average restricted to some environments within Qj\[ satisfy¬ 
ing some additional condition(s) (see Section [2.41) . We work at a fixed 
size N of the truncation domain and recall that A^(oj) is defined 
by dH). Mathematically, our approach amounts to considering condi¬ 
tional expectations of the type E[A^ | SQS], where SQS encodes that 
one, or several, of the conditions summarized in (|48p - (l49D are satisfied. 

The least we can expect from our approach is that it converges to 
the correct limit when —)• oo, namely A*, as in (j8|). 

The theorem we now state establishes this fact. In order to prove 
it, we need to make some assumptions on our setting (see the details 
below), and also to make specific the SQS conditions we use. In Theo¬ 
rem [8] below, we specifically use the SQS 1 condition, in the form (j48p . 

In order to state a result as general as possible, we therefore con¬ 


sider a condition that reads 


1 


\Qn\ 


f{Xk) = E[/(Ao)] for some 

kezd-nQN 

function /. In practice, our specific SQS 1 condition (I48p corresponds 
to the choice f{x) = x. 
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Theorem 8. Let (^fc)fcgzd be a sequence of independent and identi¬ 
cally distributed scalar random variables following a common law //. 
We assume that pL is absolutely continuous with respect to the Lebesgue 
measure on M, and that, for any k G G [—1,1] almost surely. 

We consider the stationary random field 

A{y,u}) = Co+ ^ Xkiuj)lQ+k{y)Ci{y), 


where Cq is constant and Ci is -periodic and bounded. We also 
assume that Cq + Ci{y) and Cq — C'i(y) are uniformly coercive, and 
that Cq and Ci are symmetric. 

Let f : M. M be a measurable function with compact level sets. 
We assume that f is not constant. Then we have 


E 



1 

IQa^I 


f{X,)=K[f{Xo)] 

fceQjvnZ"^ 


N^oo 


> A* 


where A^{uj) is defined by Q and A* is defined by ([7]). 
Some remarks are in order. 


(50) 


Remark 9. As is the case throughout this article, we have considered 
i/ie periodic approximation ([3]) o/ ([2]). The proof of Theorem\^ actually 
carries over to the case of Neumann or Dirichlet boundary conditions, 
or any alternate truncation problem that provides some A*'^{u;) such 
that < A*'^[uj) < j4Q^(a;) (see additional details in 

Appendix]). 

Remark 10. The assumptions regarding independence of the X^, ab¬ 
solute continuity of their common law with respect to the Lebesgue 
measure and compactness of the level sets of f are necessary for tech¬ 
nical reasons, since we need to apply a general result from J^. See 
below for details. 

The proof of Theorem [8] is based on the following result, which is 
a particular case of a more general result due to C. Bernardin and 
S. 011a (see [H Theorem B.2.2]): 

Theorem 11 (C. Bernardin and S. 011a, [3]). Consider n scalar ran¬ 
dom variables Xi, . .., Xn, that are independent and that all share 
the same probability distribution fi{x) dx on M. Consider a measurable 
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function / ; M i—>■ M, which is assumed to be not constant and to have 


compact level sets. Let /o = K[f{Xi)] = 
also a bounded and continuous function F : 


f{x)pL{x)dx. Consider 
^ !->■ M. Then 


lim E 

n^oo 


1 "" 

F{Xu...,Xk) -^f{Xi) = fo 


i=l 


= E[F{Xi,...,Xk)] 


(51) 


Note that, when n oo, the quantity — fiXAu)) almost surely 

n ^ 


i=l 


converges to /q. Theorem 1111 shows that conditioning on the manifold 
1 

— f{Xi{u)) = fo does not change the value (when n —oo) of the 


n 


2 = 1 


expectation of a function T of a finite number k of random variables. 
Note that the condition that k is independent of n can be somewhat 
relaxed. It is indeed shown in m that one can take k = o(n) in some 
cases. It is also shown there that one cannot take k = n. 


In our context, the variable X^ is the value of the field A on the 
cell Q + i. The conditioning in the left-hand side of (I5ip is identical 
to the conditioning in the left-hand side of (|50p . 

The difference between Theorem [11] and our result lies in the quan¬ 
tity of which we compute the expectation. In our case, this quantity 
is ^ 4 ^( 0 ;), which is (asymptotically when N —>■ 00 ) a function of all 
the variables Xi and not only of a finite number of them. We hence 
cannot directly use Theorem [11] The proof of our result essentially 
amounts to introducing an upper bound and a lower bound on A^(a;) 
that both read as a sum of functions that depend on a finite number 
of random variables (see e.g. (I54h below). We will then be in position 
to apply Theorem 1111 on these functions. 

Proof of Theorem\^ We fix some p G For the sake of clarity, the 
approximate homogenized matrix A^(u}) defined by (jl]) is here denoted 
Ap^ (uj), to emphasize that we have considered periodic boundary con¬ 
ditions. Since the matrix A is symmetric, we have 

P^^per{^)P = inf V € HI^^{Qn)] , 

where 

! ip + yvf A{-,u){p + Vv). 

\Qn\ jQff 
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We have considered in ([3]) periodic boundary conditions. As is well- 
known, other boundary conditions can be used, and these alternate 
approximations will be useful for the proof. 

Step 1: Upper bound. We first introduce an approximation of 
A* using a truncated corrector problem complemented with homoge¬ 
neous Dirichlet boundary conditions. We consider the problem 


-div (^A(-,a;) ^ =0 in Qat, 


w. 


N 
p,Dir 


= 0 on OQn, 


which yields an approximation of A* that we denote and which 

is defined by 




As shown in |6], we know that 


lim A 

W—>oo 


*,N 

Dir 


{u) = A* 


a.s. 


(52) 


Since A is symmetric, we have 

P^^Dhi^)P = ^ € H^{Qn)} ■ 

The matrix Apj^(a;) is always larger (in the sense of symmetric ma¬ 
trices) than Ape^(a;). Indeed, let v E Hq{Qn), and consider its Qn- 
periodic extension v. Then this function belongs to We 

hence have that 


p'^^l’^i^)P < JQNiv,u}) = Jq^{v,uj). 

Minimizing over v E HQ{Qj\f), we get that 

P^{^)P < p'^{^)P a.s. (53) 

Just as Ape^(w), the matrix Ap^(a;) depends on all the random vari¬ 
ables Xi{uj), i E Qn n But, thanks to the use of homogeneous 
Dirichlet boundary conditions, it can be bounded from above by a sum 
of matrices that depend only on a finite number of random variables. 
To show this, we proceed as follows. 
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Bn,r 




































































































jiR jiR + R 


Figure 1: The domain (here represented for iV = 11) is split into domains 
of size (here R = 2] one of them is shown in red on the figure), up to 
some boundary layer (shown in light gray). 
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For any positive integers N and R, we introduce the integer part 
M of N/R. Then Qtv can decomposed into a set of cubes of size R'^, 
up to some boundary layer 


Qn = (^ Rj + Qr) U Bn,r- 


For any j G |j| < M, consider a function vj G + Qr). 

We now define the function v on Qj\f as: 

• for any x G Rj + Qr, we set v{x) = Vj{x)‘, 

• if X G Bi\j,r, we set v{x) = 0. 

The function v belongs to Hq{Qi\[). We hence write that 


^ ■^Qn{v,uj) 


\Qr\ 

\Qn\ 


j&‘^,\j\<M 


Minimizing over the functions Vj G HQ{Rj + Qr), we hence get that 


T a*,N , N , IQrI 


^ Yj{u:) a.s. 


(54) 


\j\<M 


where 


Yj{u]) = inf {jRj+Q^{v,uj), veHl{Rj + Qr)] . 

Since A is stationary, we note that all the random variables Yj{u) 
share the same law. Moreover, we observe that lo(a;) = A]^^{oj)p, 

which is the approximation of the homogenized matrix using Dirichlet 
boundary conditions on Qr. 

We now take the conditional expectation of (|54l) . and use the fact 
that the variables Yj all share the same law: 


E 


p 


T a*,N 
^Dir 


{uj)p 


1 

IQa^I 


E 

fceQjvnZ'* 


E[/(^o)] 


<^E 

- ATd 


p 


T a*,R 
^Dir 


{uj)p 


1 

IQa^I 


E f{XQ=E[f{Xo)] 

keQNnZ'^ 


We next observe that p'^A]^l^{uj)p only depends on a finite number of 
random variables, namely only on Xk{uj) with k G Qr H Z'^. We are 
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thus in position to use Theorem [TTl which yields the limit of the above 
right-hand side when N —>• oo. Hence, for any fixed R, we have 


limsup E 

N^OO 


^ f{X,)=E[f{Xo)] 




< E 


p^A^^{uj)p 


Letting i? go to oo in the above bound and using (f52ll . we obtain that 


limsup E 

N^OO 


p^A*^^{u;)p fiXk) = nfixo)] 


< p'^ A*p. 


Using (|53p . we deduce that 


Vp G limsup 

N^oo 


(55) 


where 


Utv = E 






\Qn\ 


Y, f{Xk)=E[f{Xo)] 


fceQArnZ'^ 


(56) 


Step 2: Lower bound. We now introduce an approximation of 
A* using a truncated problem complemented with Neumann boundary 
conditions. We consider the problem 


-div (^H(-,a;) (p-h Vu;^Neu(-,w)) ) =0 in Qat, 
n^A{-,uj){p + XWp^^^^{-,uj)) =n^p on SQat, 


(57) 


which yields an approximation of A* that we denote Hj(^^(a;) and 
which is defined by 


where is dehned by 

Vp G 5*;^(w)p = J p + Vu;^Neu(-> <^)- 


(58) 


(59) 
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We refer to Remark 1121 below for some heuristic justification of (I58p - 
As recalled in [191 Appendix], we have that 

a.s. (60) 

N^oo 

and 

a.s. (61) 

In addition, we have the following variational characterization: 

= inf crGR(QAr)}, (62) 

where 

SQNi^^^) = T^[ (P + f^)'^^"H-,w)(p + cj) 

\Qn\ jQff 

and 

y{QN) = jo- G {L^{Qm)Y, divfj = 0 in Q^, n^cr = 0 on dQ^'j . 

The matrix S^^{u) (and hence the matrix A((]^(w)) depends on all 
the variables Xj(a;), i G Qn H Z'^. However, thanks to the characteri¬ 
zation (l62)) . it can be bounded from above by a sum of matrices that 
depend only on a finite number of random variables. 

To show this, we proceed as in Step 1 of the proof. For any posi¬ 
tive integers N and R, we introduce the integer part M of N/R, and 
decompose Qn into a set of cubes of size i?'^, up to some boundary 
layer Bn,r (see Figure [1]): 

Qn = ^ bijeZ'i, bj<M ^3 + Qr^ U Bn,r- 

For any j G Z'^, \j\ < M, consider a function aj G V{Rj + Qr)- We 
now define the function a on Qn as: 

• for any x G Rj + Qr, we set it(x) = crj(x); 

• if X G Bn,r, we set a{x) = 0. 

We claim that a G V{Qn)- We indeed first have that a G {L‘^{Qn))^- 
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We next consider (p G C^{Qn) and compute that 


(div cr, (f) 


-{a, V(p) 

- 

- ^ f njaj (f 

0 , 


where nj is the outward normal to the domain Rj + Qr. We hence 
have checked that a G V{Qj\f). 

We next write that 




\Qr\ 

\Qn\ 




jezd-,\j\<M 


Minimizing over the functions aj G V{Rj + Qr), we hence get that 


(63) 


where 

Zj{uj) = inf {£Rj+Q^{a,uj), a G V{Rj + Qr)} . 

Since A is stationary, we note that all the random variables Zj[u) 
share the same law. Moreover, we observe that Zq{oj) = S^^{u})p. 

We now take the conditional expectation of (|63l) . and use the fact 
that the variables Zj all share the same law: 


E 



1 

IQa^I 


f{X,) = E[fiXo)] 

keg^nzd 


j^d 

< - 

- iV<i 


P^ SNeni^)P 


1 

\Qn\ 


E 

keQisfllZd 


f{Xk) = E[/(Xo)] 


We observe that p"^S^^{(jj)p only depends on a finite number of ran¬ 
dom variables, namely only on X^ with k G Qr H We are thus 
in position to use Theorem [TTl which yields the limit of the above 
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right-hand side when N —)• oo. Hence, for any fixed R, we have 


limsup E 

N^OO 



1 

\Qn\ 


Y, fiXk)=nf{Xo)] 


< E 




Letting i? go to oo in the above bound and using (|58|l and (I60p . we 
obtain that 


limsup E 

N^OO 



1 

\Qn\ 


Y f{Xk)=nfiXo)] 




Using ()58p and (I6ip . we deduce that 


limsup E 

N^OO 


p 


, -1 




1 

\Qn\ 


E 


f{Xk) = E[/(Xo)] 




Using Jensen inequality, we infer from the above bound that 

Vp G M'^, limsup p"’" {UnY p < p'^{A*)~^p, (64) 

N^OO 

where the matrix Un is defined by ([56]). 

Step 3: Conclusion. We eventually show that (l5^ and (|MP 
imply that Un converges to A* when N ^ oo. 

From the assumptions on A, we know that there exists 0 < a_ < 
a+ < oo such that, for any N and almost surely, a_ < Hpe^(a;) < a+. 
Hence, for any N, the symmetric matrix Un satisfies a_ < Un < 
a+. We can thus extract a subsequence U(^(Ar) that converges to some 
symmetrix matrix B. Let us show that B = A*. 

Let p G We first observe that, by definition, 

limsup > lim p^U^ptj^\p = p'^Bp. 

k—>oo /c—>CXD 

We thus infer from ()55p that 

Vp G p^Bp < p"^A*p. (65) 
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We now proceed likewise with Uj^ We observe that, 

T' 1 '~r' 1 '~p 1 

limsup p UZ p > lim p UZ,^p = p B~ p. 

k^oo k^oo ‘cl, ) 

We thus infer from (I64p that 

Vp E p"^B~^p < p"^{A*)~^p. (66) 

Collecting (l65]) and ([661) . we deduce that B = A*. 

The sequence Un is bounded, and we have shown that any con¬ 
verging subsequence converges to A*. This implies that Ujsf converges 
to A* when N ^ oo, which is exactly the result (1501) . This concludes 
the proof of Theorem [S] □ 

Remark 12. In view of (j57p . we can check that 

]^\Iq =p. 

The definition (I58])-(I59]) can hence he understood as 

(^(•,w) (p +Vu;^Neu(E‘^)) ) (p +Vu;^Neu(E‘^)), 

where {■) = / • is the average on Q^. 

J Qn 

3.2 Complete analysis in some simple cases 

In this section, we aim at improving the convergence result (j50p of 
the previous section by quantifying both the statistical and systematic 
errors, in order to assess the efficiency of our approach. We are only 
able to proceed in simple situations where all the quantities are indeed 
accessible using analytic calculations. These two situations are exam¬ 
ined in Sections 13.2.11 and 13.2.21 respectively. For the sake of brevity, 
and because the proofs are not very enlightening and are not likely to 
carry over to more general cases, we do not provide the proofs of our 
claims here. We refer to |19| where they are presented in details. 

We establish below that our approach preserves the rate of decay 
of the standard Monte Carlo sampling both for the systematic and the 
statistical error (and thus, in particular, the systematic error remains, 
in rate, smaller than the statistical error). Furthermore, the prefactor 
in the statistical error is significantly reduced by our approach. 
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3.2.1 “Zero-dimensional” homogenization 


As simplest possible situation, we consider a fnnction <7 : M i—>■ M and 
the random variables (Xj)i<j<„. We assume that these random vari¬ 
ables are independent and that they are all centered Ganssian random 
variables with unit variance. We also assnme that g G C'^(M) and that 
E [|g((Aii)| -|- |gf'(Xi)|] < 00. Note that it is not snrprising to make some 
smoothness assnmptions on g as we are here after rates of convergence, 
and not only a convergence result as in Section 13.11 
We set 




1 


: X !-)• — 
n 


n 

2=1 


Assnme we want to compute E[ 5 (Xi)]. A classical Monte Carlo ap¬ 
proach would approximate this by the limit of the empirical mean 


1 -J ^ 

lim — y g{Xi{uj)). In this particular instance, the simplest version 

n—^oo n ^^ 

2=1 

1 "" 

of our variance reduction approach instead considers lim — > g{Xi(uj)) 

n^oo Ti 


2=1 


for realizations X{ijo) that satisfy ^{X{ui)) = 0 . 

In this simple case, the bias of the classical approach is actually 


identically zero: of course, E 


does not depend on n. 


n 

2=1 

The statistical error is controlled by the Central Limit Theorem and 
is asymptotically of order 


Var[5(Xi)] 


n 


Proposition 13. Under the assumptions of this section, the bias of 
the selection method is of order 1/n. More specifically, 


E 


1 

n 


2 = 1 


f{X) = 0 


-Eb(Xi)] 


-^E[g'iX,)] + 0 



(67) 

The variance of the seleetion method is reduced by a factor asymptoti¬ 
cally independent of n. More specifically, 


Var 


kT:=l9{X^) e(X)=0 


Var [iEILi 5(^0] 


Varb(Xi)] ^ 



( 68 ) 


In view of (j67p - (l68p . we observe that, at the price of introducing a 
bias of order O ( 1 /n), our approach reduces the statistical error from 
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MC _S^ (with Asqs < Amc)) and therefore, for sufficiently large 


'n y/n 
n, reduces the total error. 

The following result covers the case where we insert a non-zero 
tolerance in Algorithm [2j 

Proposition 14. Under the assumptions of this section, consider the 

selection method where we condition on the realizations such that —^ < 

y/n 

f{X{u})) < for some zq and zi > zq in R. Then, for any choice 

y/n 

of zq and zi > zq, the variance of the selection method is reduced by a 
factor asymptotically independent of n: 


Var[lEr=ig(^0 I 


Var [i Er=i 9 {Xi)] 




where C = Var 




zq<Xi< Zi 


The conditioning ze^lypn < ^(A(ti;)) < zxjypn is deliberately cho¬ 
sen in order to match the rate of the Central Limit Theorem, ft cor¬ 
responds to the selection of a fixed proportion of samples (as in Algo¬ 
rithm [3] when M is proportional to M). Note that C > 0, hence the 
variance is less reduced than when conditioning at f,{X) = 0 (which is 
the case considered in Proposition [13]) . Note also that the variance is 
reduced (with respect to the standard Monte Carlo sampling) if, and 
only if, 1 — (7 > 0. We are yet unable to conclude that this is the case 
in general. We simply note that, when zi = —zq > 0, then (7=1, 
yielding no gain. 


3.2.2 One-dimensional homogenization 

fn the one-dimensional case, the homogenization of a random field 
a : (y,uj) ^ff(Ai(a;))l(i_i+i)(y) (where g is valued, say, in [a_,a+] 
ieZ 

with a_ > 0) is a simple harmonic average, ft is readily seen that 


a*jy[{Ljj) 




with ^p{x) = l/x. 
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Formally, the problem is thus analogous to that of the previous section, 
for a certain (/? : M i—)• M instead of (/? = Id. Therefore, it is sufficient 
to prove consistency and variance reduction for quantities of the form 


Proposition 15. Consider a smooth function (/? : M i—>■ M. Under the 
assumptions of this section, the bias of the standard method and that 
of the selection method respectively are 


E 

and 




N 


N 


i=l 


-^{go) = ^^^Var[5(Xi)]+0 


iV 2 


(70) 


E 




|?W =0 


i=l 


- (ffo) 


(Var[5(Xi)] - (Eb'(Xi)])2)-^^E[Xi5'(Xi)]+o 


(71) 


with 50 = E [^(-Ei)]. 

The variance of the selection method is reduced by a factor asymp¬ 
totically independent of N: 


Var 


Ell s{x,)) 

s 

II 

o 

Var 




= 1 - 


mixiw 

Var[5(Xi)] 


+ 0 ( 1 ). (72) 


To keep things simple, we do not investigate whether a more general 
result, accounting for some tolerance in the manner our condition is 
fulfilled (in the spirit of Proposition 114)) . holds here. 

Proposition 1151 shows that the bias is unchanged in rate, while the 
prefactor for the variance is reduced. Since the variance only decays 
at the rate 1/y/N while the bias decays at the rate 1/A^, we see that 
our approach indeed reduces the total error for sufficiently large N. 

In the numerical practice (mimicking in this one-dimensional set¬ 
ting what is actually performed for higher dimensional settings - al¬ 
though it is in some sense unnecessary here), we generate several, in¬ 
dependent realizations of the A^-tuples (Xj)i<j<jv corresponding to as 
many draws of environments within the “cube” Qjsf. In the classical 
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Monte Carlo approach, we keep all such A^-tuples. In our approach, 
we only consider those that satisfy an additional criterion. 

An empirical mean (aimed at approximating A*) is then computed. 
The systematic error and the statistical error of the latter approxima¬ 
tion are precisely related to the errors estimated in (|70p - (|7ip - (|72p 
respectively. Thus a theoretical assessment of our practical approach. 


4 Numerical experiments 

We first present in this section some numerical experiments that show 
the robustness of our variance reduction approach with respect to the 
tolerance with which we enforce the SQS conditions (see Section lO) . 
We next turn to studying the performance of our approach in Sec¬ 
tion |4]2l 

We consider the test-case when A reads as in (I12p . that is 
Ar,{x,u) = Co{x,uj) -h r] xix,uj)Ci{x,uj), 
with r] = 1/2, Cq = Cl = Id, and x is of the form (j35p . that is 

X(x,w) = Xk{(jj)lQ+k{x). 
feeZ'* 

The random variables are i.i.d. and follow a Bernoulli law of pa¬ 
rameter 1/2 valued in { —1,-|-1}. The contrast (i.e. the ratio of the 
largest value of A divided by its minimum value) is equal to 3. The 
influence of the contrast on the efficiency of our approach is investi¬ 
gated at the end of Section 14.21 (see Table [T]) . We consider there much 
larger values of the contrast (however all smaller than 20 ). 

In what follows, we only consider Algorithm [3l where we take M = 
100 and A4 = 2000 (thus an acceptance ratio of 5%). 

In this setting, the SQS 1 condition as stated in (1481) is satisfied if 
and only if the numbers of cells within which Xk(uj) = 1 is equal to 
the number of cells within which Xk{oj) = — 1. It is thus possible to 
enforce ()48l) by randomly selecting \Qm\/‘2. cells within the IQa^I cells 
that are in Qtvj and setting X^ = 1 on these cells and Xk = — 1 on 
the others. 

In all our tests, we have kept the computational time fixed, or 
almost fixed, since the additional time needed by the selection step 
(namely Steps [1] and [2] of Algorithm [3]) is roughly 5% of the total 
original computational time. 
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We conclude this section with some numerical tests on a problem 
involving a more general geometry of microstructures (see Section 
On such a problem, we again obtain a significant reduction of the 
variance, at no additional computational cost. 


4.1 Robustness of the approach 


As pointed out above, the SQS 2 condition as stated in (149 p is only 
enforced in Algorithm [3] up to some tolerance. In this section, we ex¬ 
perimentally investigate how this tolerance affects the quality of the 
approximation and the efficiency of the approach. To mimick the dif¬ 
ficulty associated with the SQS 2 condition, we have also performed 
some tests where we only enforce the SQS 1 condition up to some toler¬ 
ance, and not exactly. The results of our numerical tests are displayed 
in Figures [2] through m 

Figures [2] and [3] show the sensitivity of the variance reduction ra¬ 
tio upon the first order condition (l48ll . Using Algorithm |3l we gen¬ 
erate Ai = 2000 realizations. To investigate the robustness of our 
approach, we sort these realizations with respect to the error in (148 p . 
and successively consider 20 groups of 100 realizations that less and 
less accurately satisfy (I48p . On Figure [2l the left-most circle displays 
the ratio Vsqs i/Umc between the empirical variance Vsqs i among 
the best M = 100 realizations and the reference Monte Carlo variance 
Vmc = Var 




llJ 


The second circle shows the ratio between the 
empirical variance among the next best M = 100 realizations and the 
reference Monte Carlo variance Vmc- We proceed similarly with all 
the subsequent groups of M = 100 realizations. 

On Figure [3l we display the same ratio of variances in function, 
for each group of M = 100 realizations, of the maximum error with 
which the first order condition (1481) is satisfied. Hence, the first group 
(left-most circle) corresponds to exactly satisfying the condition, the 
second group corresponds to an error between 0 and tol, the third 
group corresponds to an error between tol and 2 tol, and so on and so 
forth. 


Figures 0] and [5] show the sensitivity upon the second order condi¬ 
tion (I49p . Here, we only consider realizations that satisfy (|48p . Us¬ 
ing Algorithm [3l we again generate A4 = 2000 realizations and sort 
them according to the error in (j49p . We again successively consider 
20 groups of 100 realizations that all satisfy (1481) but that less and 
less accurately satisfy (j49p . We present the results on Figures 0] and [5] 
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following the same procedure as for Figures [2] and [3j For instance, for 

V30S 2 

the left-most circle, we plot the ratio-between the variance 

bexact SQS 1 

VsQS 2 among the M = 100 realizations that exactly satisfy the SQS 1 
condition and best satisfy the SQS 2 condition on the one hand, and, 
on the other hand, the variance 14xact SQS 1 of the realizations that 
exactly satisfy the SQS 1 condition. 

We observe that, even if the SQS conditions (I48I) - (I49I) are not ex¬ 
actly satisfied, but only with some small tolerance, we obtain a signif¬ 
icant variance reduction. We conclude that our approach is robust in 
this respect. 



Figure 2: Variance ratio Vsqs i/Vmc for the 20 groups of realizations (sorted 
according to their SQS 1 error). 

We next investigate whether enforcing the SQS 1 condition (|48p 
affects the probability distribution function of the left-hand side of the 
SQS 2 condition (j49l) . Results are shown on Figure [6l where we plot 
two histograms: 
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Figure 3: Variance ratio Vsqs i/Vmc as a function of the error in fHHj) . Results 
for only the best 7 groups (out of the 20 groups) are shown. 
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Figure 4: Variance ratio —- for the distinct groups of realizations 

Kxact SQS 1 

(sorted according to their SQS 2 error; the SQS 1 condition is exactly satis¬ 
fied). Only the 10 best groups are shown. 
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Figure 5: 
condition 


Variance ratio 


^SQS 2 
^xact SQS 1 


as a function of the error in 


fl48D is exactly satisfied). Only the 10 best groups are shown. 


(the 
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• the distribution of the criterion SQS 2 (namely, the left-hand side 
of (I49p l among all realizations. 

• the conditional distribution of the criterion SQS 2 among real¬ 
izations that exactly satisfy the SQS 1 condition (I48h . 

The two histograms are sufficiently close to each other to state that 
enforcing the SQS 1 condition does not change the distribution of the 
SQS 2 criterion. 



Figure 6: Empirical probability distribution function of the SQS 2 criterion 
(black histogram; no conditioning; red histogram: the samples exactly sat¬ 
isfy the SQS 1 criterion). Both histograms have been computed using 100 
realizations. 


4.2 Efficiency of the approach 

In this section, we investigate how the efficiency of our approach de¬ 
pends (i) on the size of the truncated domain Qpf and (ii) on the 
contrast in A. 

4.2.1 Experimental error analysis 

Figure [7] shows the set of approximations of the first entry of 

the homogenized matrix and their respective confidence intervals. We 
show three curves (along with their respective confidence intervals): 
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• The standard Monte Carlo approximation, which is defined by (uni). 
The variance is large. 

• The approximation obtained by selecting realizations that ex¬ 
actly satisfy the SQS 1 condition. The variance is much smaller, 
leading in turn to a narrower conhdence interval. 

• The approximation obtained with realizations satisfying exactly 
the SQS 1 condition and selected according to the SQS 2 condi¬ 
tion (see Algorithm [3]) . The variance is much smaller than when 
using the SQS 1 approach, even when the size of the domain Qjv 
is small. 



Figure 7; Approximations of (along with confidence intervals) as a 

function of N. Black curve: Monte Carlo method. Red curve: SQS 1 method. 
Blue curve: SQS 2 method (see text). 


Figure [8] shows a representation of the total error as a function 
of the size of Qisf. As often in a Monte Carlo approach, computing 
the total error is challenging, precisely because the reference value is, 
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by definition, in general unknown. In the specific case considered, 
namely the random checkerboard, the value of A* is actually known 
and equal to A* = {1 + r]){l — rf) Id. But for the large number of 

realizations and the large size of the truncated domains that we aim 
at considering, the total error is so small that we cannot neglect the 
contribution of the specific error due the finiteness of the meshsize used 
to solve 0. Therefore, we are bound to obtain and use as reference 
an approximate value A^ of A* corresponding to an hypothetical finite 
element approximation on the whole space. As a surrogate for this 
which is unknown in practice, we choose the empirical expectation of 
^(cj) over A4ref = 2000 random realizations exactly satisfying the 
SQS 1 condition (with a view to use a value with the lowest possible 
statistical error), and for the largest domain we can consider 

given the computing facilities we have access to, that is A^ef = 50. 

On Figure [H we display three curves: 

• The total error of the standard Monte Carlo method, defined as 


total error 


M 


M 




N 




— A* 

^ref 


m=l 


• The other two curves show the same quantity, where the M en¬ 
vironments considered now satisfy either the SQS 1 condition or 
that condition together with the SQS 2 condition. 

The results obtained using the SQS 2 approach are in general compara¬ 
ble to, and often better than, those obtained with the SQS 1 approach. 
More accurate estimates of the reference value Af^ would probably help 
in clarifying the superiority of SQS 2 over SQS 1 in terms of total accu¬ 
racy. As will now be seen, the superiority of SQS 2 in terms of variance 
(which, in some sense, is the key point for practice) is definite. 

Figure [9] shows the empirical variance of the different approxima¬ 
tions of [A*]^j^ as a function of the size of Qj\f. We again display three 
curves: 

• The standard Monte Carlo approximation defined by (|10p . 

• The approximation obtained by selecting realizations that ex¬ 
actly satisfy the SQS 1 condition. 

• The approximation obtained with realizations exactly satisfying 
the SQS 1 condition and selected according to the SQS 2 condi¬ 
tion (see Algorithmic]). 
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Figure 8: log-log plot of the total error as a function of N (natural logarithm). 
Black curve: Monte Carlo method. Red curve: SQS 1 method. Blue curve: 
SQS 2 method (see text). 
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We observe that, each time we consider an additional SQS condition, 
the empirical variance of the approximation is significantly reduced 
(even if this SQS condition is not exactly enforced; recall that we con¬ 
sider here only the 5 % best samples in terms of the SQS 2 condition, 
but that we are unable to enforce it exactly). On our test-case, en¬ 
forcing the SQS 1 condition leads to a variance 20 times smaller than 
that of the standard Monte Carlo approach, while additionally enforc¬ 
ing the SQS 2 condition leads to an additional variance reduction of a 
factor of 10. 

We also observe on Figure [9] that all variances decay as X/\Qn\, 
where 

AsQS 2 < Aexact SQS 1 < AmC- 

This corroborates in higher dimension the behaviour predicted in Sec¬ 
tion 10.21 In particular, the gain in variance does not decrease when 
the size of Qn becomes larger. 

4.2.2 Sensitivity to the contrast 

We eventually investigate how the contrast in the field A affects the 
gain in variance. Results are shown in Table [TJ We observe that the 
gain decreases when the contrast increases. Note that this is also the 
case with the antithetic variable and the control variate techniques 
that we have previously studied (see [illSlIlT]). 

However, our SQS 2 approach still yields a significant gain of a 
factor of 10 when the contrast is equal to 20. 

4.3 A more general geometry 

In order to show that the approach carries over to other settings in¬ 
volving more general geometries than the setting considered above, we 
briefly consider in the present final section a linear elasticity problem, 
for a two-phase composite material with random inclusions. The radii 
rj{u)) of the inclusions are i.i.d. random variables satisfying 



where Uj{uj) are i.i.d. variables uniformly distributed in [0,1]. The 
parameters M and m are such that the minimum (resp. maximum) 
inclusion radius is 0.125 (resp. 0.45). The inclusions centers are dis¬ 
tributed according to a Poisson point process, and we additionally 
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Figure 9: log-log plot of the variance as a function of N (natural logarithm). 
Black curve: Monte Carlo method. Red curve: SQS 1 method. Blue curve: 
SQS 2 method. 
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Contrast 

Vmc 

^xact SQS 1 

fosQS 2 

Vmc 

Vmc 

^xact SQS 1 

fosQS 2 

1.22 

0.0000273 

5.801e-08 

6.858e-10 

470 

39821 

1.50 

0.0001097 

0.0000009 

1.585e-08 

118 

6921 

1.86 

0.0002488 

0.0000047 

0.0000001 

52.6 

1996 

2.33 

0.0004478 

0.0000151 

0.0000006 

29.5 

720 

3.00 

0.0007118 

0.0000379 

0.0000024 

18.8 

296 

4.00 

0.0010496 

0.0000814 

0.0000080 

12.8 

131 

5.67 

0.0014769 

0.0001600 

0.0000244 

9.23 

60.5 

9.00 

0.0020289 

0.0003021 

0.0000739 

6.71 

27.4 

19.0 

0.0028330 

0.0006061 

0.0002554 

4.67 

11.1 


Table 1: For various values of the contrast, we show the Monte Carlo vari¬ 
ance (column ^ 2 ), the variance of the SQS 1 method (column ^3) and the 
variance of the SQS 2 method (column 7 ^ 4 ). We next show the variance ratio 
hMc/Kxact SQS 1 for fhe SQS 1 approach (column 7^5) and the variance ratio 
VMc/fosQS 2 for the SQS 2 approach (column 7 (^ 6 ). The size of Qn is hxed at 
N = 20. 


impose that inclusions do not overlap (see Figure EOD. We consider 
the truncated domain Qj^ = [0,5]^. Each microstructure contains 25 
inclusions. If some part of an inclusion falls outside of it is repro¬ 
duced on the other side of Qn by periodicity. 

The inclusions (resp. the background) are modeled by a isotropic 
linear elasticity tensor, with a uniform Poisson ratio ly = 0.3 and a 
Young modulus E = 7 (resp. E = 1). 

For that problem, in the spirit of the SQS 1 approach described 

above, we select the microstructures such that the volumic fraction 

0n{oj) of inclusions within the domain Qn-, for each random realization 

considered, agrees as accurately as possible with its asymtotic value 

r= lim 07 v(a;). 

N—^OO 

The results are shown in Table [2j We examine two entries of the 
homogenized elasticity tensor, and two methods: 

• The classical Monte Carlo approach, for which we generate M = 
100 i.i.d. microstructures. 

• The SQS approach, in which we generate Ai = 2000 i.i.d. mi¬ 
crostructures, and next consider the M = 100 microstructures 
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Figure 10; Two realizations of the microstructure geometry. 


for which \6i^{uj) — 9*\ is the smallest. 

In both cases, we solve M correctors problems and we compute the 
empirical expectation and variance of the homogenized elasticity ten¬ 
sor. We observe that our approach provides a variance reduction of a 
factor close to 20, while the bias is essentially constant. 



[y4^(a;)]iiii 

[y4^(a;)]iiii 

[y4^(a;)]ii22 

[y4^(a;)]ii22 


Exp. 

Var. 

Exp. 

Var. 

MC approach 

2.522 

0.0136 

1.016 

0.00184 

SQS approach 

2.519 

0.000532 

1.014 

0.000101 


Table 2: Empirical expectation (Exp.) and variance (Var.) for two entries of 
the homogenized elasticity tensor, computed with the Monte Carlo approach 
(MC) or our approach (SQS). 
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A 


Proof of Lemma 


a 


We follow the arguments of the proof of [5l Lemma 3.2]. 

The existence and uniqueness (up to the addition of a constant) of 
(j)i solution to (1381) is established in [5l Lemma 3.1]. We next point out 
that fl39p admits a unique (up to the addition of a constant) solution 
in iLpgj.((5). It is a simple consequence of the Lax-Milgram lemma. 

We now prove that the sum in (I37p is a convergent series in {Q x 
D). For this purpose, we compute the norm of the remainder of the 
series, using the notation X^ioj) = X^ioj) — E[Afc]: 

2 


^ XkV4>i{--k) 

|fc|>Af-|-l 


L2(Qxn) 


= ^ ^ / Vct>i{y-k)-V^{y-l)dy 

|fc|>Af-|-l |£|>Ar-|-l •'Q 


< E E |Cov(Afc,A,)| \\V4>i\\l^^q-i) 

|fc|>Ar-|-l |£|>Af-|-l 


< E E |Cov(Afc,A,)| 

|fc|>7V-|-l |£|>Ar+l 


where we have used at the last line the discrete Gauchy-Schwarz in¬ 
equality between the sequences |Cov(Afc, and 

|Cov(Afc, \\X We next write, using the stationarity 

of Afc and ([36l) . that 


< 


< 


^ XkV4>i{--k) 

L2(Qxr2) 

E l|V0i|li2(Q_,) |Cov(Afc,A, 

|fc|>Af-|-l |£|>Af-|-l 
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The above right-hand side converges to 0 as —>■ oo since V(/>i E 

Hence, the right-hand side of (I37p defines a function T E {L‘^{Q x 
As diTj = djTi, there exists a function ui such that 

vSi = r = E[Xo]Vur+ ^ -E[Xk])vM- - k). 

As uT is Z'^-periodic, we infer from the above equality that 


Vtti is stationary and 



E(Vui) = 0. 


Next, we compute 


(73) 


CoXui = E[Xo]CoVur+ [Xkiu] - E[Xk])CoXM- " k). 


Taking the divergence of this equation and using (1321) and (|34p . we 
thus find that, in the distribution sense, 


— div [CoVui] 

= - [Xk{u) - E[Afc]) div [CoVM- - k)] - E[Ao] div [CqVW] 

fceZ'* 

= ^ (Afc(w) - E[Afc]) div [iQ+kCip] + E[Xo] div [Cip] 

= div[x(-,w)C'ip]. (74) 

CollGcting (jT3p Siiici f|T4:P . wg sgg tlicit soIvgs m- As the solution 
to this equation is unique up to the addition of a (possibly random) 
constant C{oj), we obtain that tti = tti -|- C{(jo), hence proving (|37p . 
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